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ABSTRACT 


A systematic  method  of  constructing  formulae,  together  with  error 
terms,  is  given  for  use  in  interpolation,  extrapolation,  differentiation, 
and  integration.  Both  closed  and  open  type  formulae  are  developed,  using 
ordinates,  and  based  on  Lagrange’s  interpolation  formulae  for  equal  inter- 
vals.  The  procedure  was  suggested  by  Professor  H.  H.  Aiken  several  years 
ago  at  Harvard  University.  The  whole  procedure  may  be  extended  easily  to 
obtain  cubature  formulae  and  formulae  used  for  surface  fitting. 
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LAGRANGIAN  INTERPOLATION 


The  purpose  of  this  paper  is  to  apply  Lagrange's  interpolation 
formulae  1,  2,  3,  4 to  derive  polynomial  approximating  formulae  (using 
ordinates)  for  interpolation,  differentiation,  and  integration.  Many 
of  these  formulae  appear  in  the  literature®  in  terms  of  differences.  Of 
course  these  differences  may  be  expressed  in  terms  of  ordinates  and  in  a 
sense  the  formulae  which  appear  here  are  a recapitulation  of  elementary 
formulae  fcund  in  Numerical  Analysis.  However,  the  simple  derivations 
given  here  admit  easy  extensions  to  new  formulae  of  great  accuracy  (ex- 
clusive of  round-off  error)  to  be  used  in  connection  with  large  scale 
high  speed  digital  calculators.  The  number  of  ordinates  employed  in  the 
formulae  are  indicated  in  the  tables  and  include  as  many  as  ten. 

LA(3RANGE'S  INTERPOUTION  FORMULA 

The  following  assumptions  are  made  and  the  notation  is  more  or  less 
standard. 
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(a)  X = Xo  + uh,  y,  = y(xo  + sh)  = y(x,). 

(b)  s,  a,  b are  integers. 

(c)  a < u < b,  where  u may  be  integral  or  not. 

(d)  n + 1 = b - a + 1 is  the  number  of  points  through  which  the 
approximating  polynomial  (1)  passes. 

(e)  £ lies  between  the  smallest  and  the  largest  of  the  set  x, 

where  s - a = 0,  2,...,  n. 

(f)  (£)  is  familiar  notation  for  the  (n  + l)st  derivative 
evaluated  at  x = £. 

Lagrange's  formula  may  be  written 

b 

(1)  y(x)  = ^ C,.,(x)  y(x)  + R(x)  where 
s=a 


C (v)  - ■ - ' 

(x-x.)<^'(x.)  (u-s)<^'(s) 


, and 


(2)  R(x)  = 


<^(u)  y(n^i)  (g)  ^ 
(n  + 1).' 


Also  by  definition 


<^(x)  * n (x  + ah  - X,)  = h"'^^  n (u  - s)  = h"'*'^ 

s=0 


s=a 


Proof  of  (1). 

X ) 

Assume  y = \ v.  + R(x). 

Z_  x-x, 
s=a 


s (an  integer)  we  haye 


Then  for  any  particular 
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1-  ^ tl 

y.  = V-  lim = (X.).  Hence 

“ " x-3^  X - x^  * “ 

ys 

V.  = — TT — r which  leads  to  (1). 

(x,) 

Proof  ® of  (2). 

Write  f(z)  = R(z)  - R(x)  77—  • It  follows  that  f = 0 at  z = x,  x 

<^(x)  « 

for  s - a = 0,  1,  2,  3,...,  n.  Hence  by  Rolle's  theorem  ’’  ) = 0. 

Also  it  is  apparent  that  ^^“‘‘‘^^(z)  = (n  + 1).'  Consequently 

0 = R(“^i)  (£)  - R(x)  so  that  R(x)  = <^(x)  R^^+l)  (£) 

<^(x)  — 

(n  + 1).' 

Now  by  (1) 

R(n+i)  (g)  , y(n+i>(g)  which  proves  (2). 

Ifb<uoru<a  then  (1)  may  be  used  to  extrapolate  y if  we  accept 
the  error  produced  by  taking  t in  (2)  to  lie  in  the  interval  between  the 
least  and  greatest  of  the  set  x,  x,  and  again  where  s - a = 0,  i,  2,...,  n 


DIFFERENTIATION  FORMULAE 


By  differentiating  (1)  and  substituting  x = x^  where  t is  a par- 
ticular s we  get 


(3) 


C'.^y^  + R'(t)  where 
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(4)  R'(t) 


h“<^'(t)  <£) 

(n  + 1): 


In  order  to  prove  (4)  consider  (2)  which  may  be  written 


y(n+l)(g) 

R(x)  = ^6(x)  y(x,  X.)  since  is  a function  of  the  set  x,  x 

® (n  + 1)J  ® 

where  s - a = 0,  1,  2,...,  n.  Hence  R'  = 4'y  + 4y‘  and  since 

<^(x  ) = 0 and  <^'(x)  = h“'*'^  ^ ^ = h“  <^'(u)  we  have  (4). 

^ du  dx 

In  order  to  simplify  the  remainder  terms  the  derivative  formulae 
are  tabulated  at  the  given  ordinates  only.  Moreover  not  all  these  are 
tabulated  since  some  may  be  implied.  For  example  the  29  first  deriva- 
tive formulae  imply  25  more.  Indeed  for  s / t we  have 


c;  (t)  * ~ lira 


<^(u) 


•77-  , and 


K u-t  <^1 

c;..  („-  t)  . Im  ^ . lim 


U-71-S  u-n-t  (a-n+t)  (u-n+s) 


= lira 


u-n+s 


. lim 


^(n-u) 


u-Ti-s  <^(n-u)  u n-t  (u-n+t)  (u-n+s) 


= - lira 


u-s 


lira 


^(u) 


K 


<^(u)  ' U-t  (u-t)  (u-s)  " t-s  * r ’ 

Hence  C^(t)  = — (n  — t).  The  diagonal  elelnents  obey  this  relation 


also  since 
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n 


z 

s=0 


(t) 


n 

L 


C^.(n  - t) 


0. 


In  order  to  obtain  higher  order  differentiation  formulae  we  may 
differentiate  y*  " ^ ys  » the  remainder  for  a moment. 

Here  y'  and  y^  are  vectors  while  is  a square  matrix  of  order  n + 1. 
Hence  in  general 

(5)  y<"^  = C'"  y . 

The  remainder  is  obtained  by  differentiating  R(x)  = <^(x)  y Cx,  x,) 

as  before.  We  have  in  general,  using  symbolic  notation, 

(6)  R^"^(x)  = [^(x)  + y(x,  x,)]^"^ 

Only  the  first  term  of  R^"^(x)  is  retained  i.e.  R^"^(x)  = <^^“^(x)  y(x,  x,) 
unless  this  term  vanishes.  If  <^^“^(x)  = 0 then  obviously  the  order  of 
R^"^(x)  is  multiplied  by  h and  is  so  indicated. 

The  implied  differentiation  formulae  are  obtained  by  using  the 
relation  Ci">  (t)  = C<“>(n-t). 

s n * s 

In  the  tabulation  below  r^“^  = - — 7— • 

(n+l)J 
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QUADRATURE  FORMULAE 

We  may  obtain  quadrature  formulae  by  integrating  (1)  to  get 


(7) 


b+1 


b-k 


ydx  “ ^ ^ ^s-a  ^a  where 


s=a 


},n+2  „(n+l)  tp)  f 1 

(8)  R*  < f - <^{U  + n)  du 

(n  + 1)J  . 


(unless  the  last  integral  vanishes)  and  where 

n+1  <^(u)  du 

+ nj  du  = — I 

At  J T 

K n-k 


c;  . I C.(u  t n)  du  - ^ I 


u - s 


In  order  to  prove  (8)  consider  again 

R(x)  = ^(x)  y(x,  x^)  from  which 
rl 

R*  = h’*'*'*  <^(u  + n)  y(x,  Xg)  du 

\*-k 

n 

Now  <^(u  + n)  = n (u  + s)  which  is  raonotonic  increasing  for  u > 
s=0 

By  the  Mean  Value  Theorem  we  may  therefore  write 

(^1 ) 

<^(u  + n)  y(x,  X ) du  = ■■  , — —7-  <^(u  + n)  du  where 

Jo  (n  + !)•  Jo 

0 < 5i  < 1.  Again  by  Rolle's  Theorem  since  <^(u  + n)  is  continuous  and 
vanishes  at  u = 0,  -1,  there  is  a point  -1  < t?  < 0 such  that  <^'(u  + n) 
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Consequently  in  the  interval  17  < u < 0,  <^(u  + n)  is  monotonic  increasing 
so  that  we  may  write 


4in  + n)  y(x,  du  = 


(n  + 1).’ 


<^(u  + n)  du  where 


7)  < £2  '^0.  If  y(“'*'i)  (£)  is  the  larger  of  y<“'*’i)  (£1),  y^“'’’^U£2)  then 
we  may  write 


+ n)  y(x, 


X, ) du  < 


(n  + )1.» 


+ n)  du  where 


7J  < £ < 1. 

Continuing  for  the  interval  -1  < u < 17  we  find  that  <^(u  + n)  is 
monotonic  decreasing  so  that  we  may  write 


■77  y(n+l)(£  ) ,77 

<^(u  + n)  y(x,x.)  du  * . . , ^(u  + n)  du  where 

-1  • (n  * 1);  J_^ 


-!<£,<  77.  Therefore  we  may  write 


+ n)  y(x,  x^)  du 


y<“"-l>  (£2)  I + y(n+l)  (gj 

77  • 


y(«+l)  (£) 

■ (n+1)' 


u + n)  du 


where  -1  < £ < 1 and  y^“‘*'t)  (£)  is  the  largest  of  C£i),  (£2), 

y(n+i)  (gg).  Since  the  polynomial  vanishes  at  -u  = 0,  1,  2,...,  n we 
find  (8)  by  continuing  in  this  way. 
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The  condition  <^(u  + n)  du  =7^  0 cannot  be  relaxed.  However  because 
-k 

of  this  condition  no  formula  will  be  lost.  For  consider 


■#.-1,  -1  c;.,. 

<^(u) 

du 

<^(u-l)  du 

n-l  = J 

n-k-1 

(u-s+1 ) 

(u-n) 

•J  (u-s)  (u-n-1) 

n-k 

l^n+1 

“^n-k 

^(u)  du 
u(u-s) 

. Hence 

we  have 

c:  - s 

n-l  ^^-1,  n-l 

/-n+1 

n-k 

[ _ 
u-s 

- = *'“>  ] du 

u(u-s) 

|-n  + l 

n-k 

^(u)  du 
u 

<ct 

and  since  s ^ we  have 

1 ' ■ 


(9)  c:  = ^ cs  . c:.i, 


The  C*_j^  refers  to  that  value  immediately  above  C*.i  in  the  tables 

below.  Of  course 
rl 

Cj  =1  (u  + n)  du  where  the  integral  is  on  the  line 

-k 


immediately  above  Cj  in  the  tables.  If  this  integral  vanishes  we  have 


^8  ■ ^8-1,  n-i  s y 0,  and  therefore  the  formula  is  repeated. 
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It  is  easy  to  see  that  this  integral  cannot  vanish  twice  in  succession. 

By  taking  1 = 1 we  obtain  open  quadrature  formulae  and  by  taking 
1 < 0 we  obtain  closed  quadrature  formulae.  Because  of  symmetry  the 
implied  quadrature  formulae  double  to  number  of  formulae  listed  below. 

For  example  if 'C*  is  replaced  by  C*.^  in  a formula,  for  extrapolating 

ahead  we  shall  replace  it  by  one  which  extrapolates  backward.  The 
remainders  may  be  different  because  the  intervals  are  different  for  £. 
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APPLICATION  OF  FORMULAE 


Desirable  closed  formulae  are  easily  seen  to  be  No's  64,  78,  88, 
and  94.  Some  caution  should  be  used  when  the  8-strip  formula  is  employed 
for  since  some  of  its  coefficients  are  negative  the  round-off  error  may 
appear  excessive.  Although  some  of  the  coefficients  are  negative,  nevei^ 
the-less,  open  formulae  No's  27,  40,  49,  and  54  show  considerable  appeal. 

For  example  No.  40  was  employed  in  A.F.  Problem  No.  19  at  Harvard  Con- 
putation  Laboratory. 

Very  little  can  be  said  in  favor  of  increasing  the  error  in  order 
to  simplify  the  coefficients  when  the  numerical  operations  are  performed 
on  a large  scale  digital  electronic  calculator.  Consequently  Weddle's 
Rule,  Hardy's  Formula,  Shovelton's  Rule,  etc.  will  not  be  introduced  here. 

Often  it  is  necessary  to  change  the  mesh  interval  when  integrating 
a system  of  differential  equations.  Formula  (1)  may  be  used  to  generate 
a set  of  formulae  to  cut  such  an  interval.  Obviously  no  special  formulae 
are  necessary  to  double  an  interval. 

It  is  interesting  to  note  that  by  decreasing  the  mesh  size  indefinitely 
and  at  the  same  time  keeping  all  the  ordinates  (an  infinite  number  in  a 
finite  interval)  we  are  led  to  the  cardinal  interpolation  function  mentioned 
in  AFTR  No.  65581  by  P.  W.  Bubb. 
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